!#-----------------------------------------------------------------------------#
!# This file was copied from an early version from GOTM and modified to remove #
!# the explicit array boundary information in routine arguments, then          #
!# modularised and support added for single precision compilation              #
!#                                                                             #
!#  Original author(s): Hans Burchard, Karsten Bolding, Jorn Bruggeman         #
!#  Copyright by the GOTM-team under the GNU Public License - www.gnu.org      #
!#                                                                             #
!#                   ------------------------------                            #
!#                                                                             #
!#  These mods :                                                               #
!#     AquaticEcoDynamics (AED) Group                                          #
!#     School of Agriculture and Environment                                   #
!#     The University of Western Australia                                     #
!#                                                                             #
!#     http://aquatic.science.uwa.edu.au/                                      #
!#                                                                             #
!# Copyright 2013 - 2022 -  The University of Western Australia                #
!#                                                                             #
!#  This file is part of GLM (General Lake Model)                              #
!#                                                                             #
!#  GLM is free software: you can redistribute it and/or modify                #
!#  it under the terms of the GNU General Public License as published by       #
!#  the Free Software Foundation, either version 3 of the License, or          #
!#  (at your option) any later version.                                        #
!#                                                                             #
!#  GLM is distributed in the hope that it will be useful,                     #
!#  but WITHOUT ANY WARRANTY; without even the implied warranty of             #
!#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              #
!#  GNU General Public License for more details.                               #
!#                                                                             #
!#  You should have received a copy of the GNU General Public License          #
!#  along with this program.  If not, see <http://www.gnu.org/licenses/>.      #
!#                                                                             #
!#-----------------------------------------------------------------------------#

#include "glm.h"

# define REALTYPE AED_REAL

!#define _ARRSZ_ 1:numc,1:nlev
!#define _ARRIDX_ ci,i

#define _ARRSZ_ 1:nlev,1:numc
#define _ARRIDX_ i,ci


#define _POINT_NEG_NINE_ 1.d-9

MODULE ode_solvers

   USE ISO_C_BINDING

   IMPLICIT NONE

CONTAINS

#if 1
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: General ODE solver \label{sec:ode-solver}
!
! !INTERFACE:
   subroutine ode_solver(solver,numc,nlev,dt,cc,right_hand_side_rhs,right_hand_side_ppdd)
!
! !DESCRIPTION:
! Here, 10 different numerical solvers for the right hand sides of the
! biogeochemical models are implemented for computing the ordinary
! differential equations (ODEs) which are calculated as the second
! step of the operational split method for the complete biogeochemical
! models. The remaining ODE is
! \begin{equation}\label{ODESystem}
!   \partial_t c_i
! =   P_i(\vec{c}) -D_i(\vec{c}), \;\; i = 1,\ldots,I,
! \end{equation}
! with $c_i$ denoting the concentrations of state variables.
! The right hand side denotes the reaction terms,
! which are composed of contributions
! $d_{i,j}(\vec{c})$, which represent reactive fluxes from
! $c_i$ to $c_j$, and in turn, $p_{i,j}(\vec{c})$ are reactive fluxes from
! $c_j$ received by $c_i$, see equation (\ref{eq:am:a}).
!
! These methods are:
!
! \begin{enumerate}
! \item First-order explicit (not unconditionally positive)
! \item Second order explicit Runge-Kutta (not unconditionally positive)
! \item Fourth-order explicit Runge-Kutta (not unconditionally positive)
! \item First-order Patankar (not conservative)
! \item Second-order Patankar-Runge-Kutta (not conservative)
! \item Fourth-order Patankar-Runge-Kutta (does not work, not conservative)
! \item First-order Modified Patankar (conservative and positive)
! \item Second-order Modified Patankar-Runge-Kutta (conservative and positive)
! \item Fourth-order Modified Patankar-Runge-Kutta
!       (does not work, conservative and positive)
! \item First-order Extended Modified Patankar
!       (stoichiometrically conservative and positive)
! \item Second-order Extended Modified Patankar-Runge-Kutta
!       (stoichiometrically conservative and positive)
! \end{enumerate}
!
! The schemes 1 - 5 and 7 - 8 have been described in detail by
! \cite{Burchardetal2003b}. Later, \cite{Bruggemanetal2005} could
! show that the Modified Patankar schemes 7 - 8 are only conservative
! for one limiting nutrient and therefore they developed the
! Extended Modified Patankar (EMP) schemes 10 and 11 which are also
! stoichiometrically conservative. Patankar and Modified Patankar
! schemes of fourth order have not yet been developed, such that
! choices 6 and 9 do not work yet.
!
! The call to {\tt ode\_solver()} requires a little explanation. The
! first argument {\tt solver} is an INTEGER and specifies which solver
! to use. The arguments {\tt numc} and {\tt nlev} gives the dimensions
! of the data structure {\tt cc} i.e. {\tt cc(1:numc,0:nlev)}.
! {\tt dt} is simply the time step. The last argument is the most
! complicated. {\tt right\_hand\_side} is a subroutine with a fixed
! argument list. The subroutine evaluates the right hand side of the ODE
! and may be called more than once during one time-step - for higher order
! schemes. For an example of a correct {\tt right\_hand\_side} have a look
! at e.g. {\tt do\_bio\_npzd()}
!
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   INTEGER, intent(in)                 :: solver,nlev,numc
   REALTYPE, intent(in)                :: dt
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side_ppdd(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)

         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end
   end interface

   interface
      subroutine right_hand_side_rhs(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------
!BOC
   select case (solver)
      case (1)
         call euler_forward(dt,numc,nlev,cc,right_hand_side_rhs)
      case (2)
         call runge_kutta_2(dt,numc,nlev,cc,right_hand_side_rhs)
      case (3)
         call runge_kutta_4(dt,numc,nlev,cc,right_hand_side_rhs)
      case (4)
         call patankar(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (5)
         call patankar_runge_kutta_2(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (6)
         call patankar_runge_kutta_4(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (7)
         call modified_patankar(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (8)
         call modified_patankar_2(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (9)
         call modified_patankar_4(dt,numc,nlev,cc,right_hand_side_ppdd)
      case (10)
         call emp_1(dt,numc,nlev,cc,right_hand_side_rhs)
      case (11)
         call emp_2(dt,numc,nlev,cc,right_hand_side_rhs)
      case default
         stop "glm_fabm: no valid solver method specified in glm.nml !"
   end select

    ! STDERR 'ode_solver return'
   return
   end subroutine ode_solver
!EOC
#else
SUBROUTINE init_ode_solver(ode_method, solver)
  INTEGER,INTENT(in) :: ode_method
  PROCEDURE(),pointer,INTENT(out) :: solver
!BEGIN
   !# Report type of solver
   print *, "Using Eulerian solver"
   SELECT CASE (ode_method)
      CASE (1)
         print *, 'Using euler_forward()'
         solver => euler_forward ! (dt,numc,nlev,cc,right_hand_side_rhs)
      CASE (2)
         print *, 'Using runge_kutta_2()'
         solver => runge_kutta_2 ! (dt,numc,nlev,cc,right_hand_side_rhs)
      CASE (3)
         print *, 'Using runge_kutta_4()'
         solver => runge_kutta_4 ! (dt,numc,nlev,cc,right_hand_side_rhs)
      CASE (4)
         print *, 'Using patankar()'
         solver => patankar ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (5)
         print *, 'Using patankar_runge_kutta_2()'
         solver => patankar_runge_kutta_2 ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (6)
         print *, 'Using patankar_runge_kutta_4()'
         solver => patankar_runge_kutta_4 ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (7)
         print *, 'Using modified_patankar()'
         solver => modified_patankar ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (8)
         print *, 'Using modified_patankar_2()'
         solver => modified_patankar_2 ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (9)
         print *, 'Using modified_patankar_4()'
         solver => modified_patankar_4 ! (dt,numc,nlev,cc,right_hand_side_ppdd)
      CASE (10)
         print *, 'Using emp_1()'
         solver => emp_1 ! (dt,numc,nlev,cc,right_hand_side_rhs)
      CASE (11)
         print *, 'Using emp_2()'
         solver => emp_2 ! (dt,numc,nlev,cc,right_hand_side_rhs)
      CASE DEFAULT
         STOP 'init_glm_fabm: no valid ode_method specified in glm.nml!'
   END SELECT

END SUBROUTINE init_ode_solver
#endif

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: First-order Euler-forward scheme
!
! !INTERFACE:
   subroutine euler_forward(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the first-order Euler-forward (E1) scheme is coded, with one
! evaluation of the right-hand sides per time step:
! \begin{equation}\label{eq:am:euler}
! \begin{array}{rcl}
! c_i^{n+1} &=&
! \displaystyle
!  c_i^n  + \Delta t \left\{P_i\left(\underline{c}^n\right)
! - D_i\left(\underline{c}^n\right) \right\}.
! \end{array}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   INTEGER, intent(in)                 :: numc,nlev
   REALTYPE, intent(in)                :: dt
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  REALTYPE :: rhs(_ARRSZ_)
  INTEGER  :: i,ci
!EOP
!-----------------------------------------------------------------------
!BOC
   call right_hand_side(.true.,numc,nlev,cc,rhs)

   do i=1,nlev
      do ci=1,numc
         cc(_ARRIDX_)=cc(_ARRIDX_)+dt*rhs(_ARRIDX_)
      end do
   end do

   return
   end subroutine euler_forward
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Second-order Runge-Kutta scheme
!
! !INTERFACE:
   subroutine runge_kutta_2(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the second-order Runge-Kutta (RK2) scheme is coded, with two
! evaluations of the right hand side per time step:
! \begin{equation}\label{eq:am:RK}
! \left.
! \begin{array}{rcl}
! c_i^{(1)} &=&
! \displaystyle
! c_i^n  + \Delta t \left\{P_i\left(\underline{c}^n\right)
! - D_i\left(\underline{c}^n\right) \right), \\ \\
!  c_i^{n+1} &=&
! \displaystyle
!  c_i^n  + \dfrac{\Delta t}{2}
! \left\{P_i\left(\underline{c}^n\right) + P_i\left(\underline{c}^{(1)}\right)
!  - D_i\left(\underline{c}^n\right) - D_i\left(\underline{c}^{(1)}\right)
! \right\}.
! \end{array}
! \right\}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  REALTYPE :: rhs(_ARRSZ_),rhs1(_ARRSZ_)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,ci
!EOP
!-----------------------------------------------------------------------
!BOC
   call right_hand_side(.true.,numc,nlev,cc,rhs)

   do i=1,nlev
      do ci=1,numc
         cc1(_ARRIDX_)=cc(_ARRIDX_)+dt*rhs(_ARRIDX_)
      end do
   end do

   call right_hand_side(.false.,numc,nlev,cc1,rhs1)

   do i=1,nlev
      do ci=1,numc
         cc(_ARRIDX_)=cc(_ARRIDX_)+dt*0.5*(rhs(_ARRIDX_)+rhs1(_ARRIDX_))
      end do
   end do

   return
   end subroutine runge_kutta_2
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Fourth-order Runge-Kutta scheme
!
! !INTERFACE:
   subroutine runge_kutta_4(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the fourth-order Runge-Kutta (RK4) scheme is coded,
! with four evaluations
! of the right hand sides per time step:
! \begin{equation}\label{eq2}
! \left.
! \begin{array}{rcl}
! c_i^{(1)} &=&
! \displaystyle
! c_i^n+\Delta t \left\{P_i\left(\underline c^n\right)
! -D_i\left(\underline c^n\right)\right\} \\ \\
! c_i^{(2)} &=&
! \displaystyle
! c_i^n+\Delta t \left\{P_i\left(\frac12\left(\underline c^n+\underline
! c^{(1)}\right)\right)-D_i\left(\frac12\left(\underline c^n+\underline
! c^{(1)}\right)\right)\right\} \\ \\
! c_i^{(3)} &=&
! \displaystyle
! c_i^n+\Delta t \left\{P_i\left(\frac12\left(\underline c^n+\underline
! c^{(2)}\right)\right)-D_i\left(\frac12\left(\underline c^n+\underline
! c^{(2)}\right)\right)\right\} \\ \\
! c_i^{(4)} &=&
! \displaystyle
! c_i^n+\Delta t \left\{P_i\left(\underline c^{(3)}\right)-D_i\left(\underline
! c^{(3)}\right)\right\} \\ \\
! c_i^{n+1} &=&
! \displaystyle
!  \frac16 \left\{c_i^{(1)}+2c_i^{(2)}+2c_i^{(3)}+c_i^{(4)}   \right\}.
! \end{array}
! \right\}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: rhs(_ARRSZ_),rhs1(_ARRSZ_)
  REALTYPE :: rhs2(_ARRSZ_),rhs3(_ARRSZ_)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,ci
!EOP
!-----------------------------------------------------------------------
!BOC
   first=.true.
   call right_hand_side(first,numc,nlev,cc,rhs)
   first=.false.

   do i=1,nlev
      do ci=1,numc
         cc1(_ARRIDX_)=cc(_ARRIDX_)+dt*rhs(_ARRIDX_)
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,rhs1)

   do i=1,nlev
      do ci=1,numc
         cc1(_ARRIDX_)=cc(_ARRIDX_)+dt*rhs1(_ARRIDX_)
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,rhs2)

   do i=1,nlev
      do ci=1,numc
         cc1(_ARRIDX_)=cc(_ARRIDX_)+dt*rhs2(_ARRIDX_)
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,rhs3)

   do i=1,nlev
      do ci=1,numc
         cc(_ARRIDX_)=cc(_ARRIDX_)+dt*1./3.*(0.5*rhs(_ARRIDX_)+rhs1(_ARRIDX_)+rhs2(_ARRIDX_)+0.5*rhs3(_ARRIDX_))
      end do
   end do

   return
   end subroutine runge_kutta_4
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: First-order Patankar scheme
!
! !INTERFACE:
   subroutine patankar(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the first-order Patankar-Euler scheme (PE1) scheme is coded,
! with one evaluation of the right hand sides per time step:
! \begin{equation}\label{eq:am:patankar}
! \begin{array}{rcl}
!   c_i^{n+1} &=&
! \displaystyle
!  c_i^n  + \Delta t \left\{P_i\left(\underline{c}^n\right)
!                                       - D_i\left(\underline{c}^n\right)
!                                         \frac{c_i^{n+1}}{c_i^n} \right\}.
! \end{array}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: ppsum,ddsum
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_
   dd=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do i=1,nlev
      do ci=1,numc
         ppsum=_ZERO_
         ddsum=_ZERO_
         do j=1,numc
            ppsum=ppsum+pp(i,j,ci)
            ddsum=ddsum+dd(i,j,ci)
         end do
         cc(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum)/(1.+dt*ddsum/cc(_ARRIDX_))
      end do
   end do

   return
   end subroutine patankar
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Second-order Patankar-Runge-Kutta scheme
!
! !INTERFACE:
   subroutine patankar_runge_kutta_2(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the second-order Patankar-Runge-Kutta (PRK2) scheme is coded,
! with two evaluations of the right hand sides per time step:
!
! \begin{equation}\label{eq:am:PRK}
!   \left.
!   \begin{array}{rcl}
!     c_i^{(1)} &=&
! \displaystyle
!  c_i^n  + \Delta t
!                   \left\{P_i\left(\underline{c}^n\right)
!                       - D_i\left(\underline{c}^n\right)
!                         \dfrac{c_i^{(1)}}{c_i^n}\right\},
!                   \\ \\
!     c_i^{n+1} &=&
! \displaystyle
!  c_i^n  + \dfrac{\Delta t}{2}
!                   \left\{P_i\left(\underline{c}^n\right)
!                         + P_i\left(\underline{c}^{(1)}\right)
!                         - \left( D_i\left(\underline{c}^n\right)
!                         + D_i\left(\underline{c}^{(1)}\right)\right)
!                         \dfrac{c_i^{n+1}}{c_i^{(1)}}
!                   \right\}.
!   \end{array}
!   \right\}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: ppsum(_ARRSZ_),ddsum(_ARRSZ_)
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_
   dd=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do ci=1,nlev
      do i=1,numc
         ppsum(_ARRIDX_)=_ZERO_
         ddsum(_ARRIDX_)=_ZERO_
         do j=1,numc
            ppsum(_ARRIDX_)=ppsum(_ARRIDX_)+pp(i,j,ci)
            ddsum(_ARRIDX_)=ddsum(_ARRIDX_)+dd(i,j,ci)
         end do
         cc1(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum(_ARRIDX_))/(1.+dt*ddsum(_ARRIDX_)/cc(_ARRIDX_))
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,pp,dd)

   do ci=1,nlev
      do i=1,numc
         do j=1,numc
            ppsum(_ARRIDX_)=ppsum(_ARRIDX_)+pp(i,j,ci)
            ddsum(_ARRIDX_)=ddsum(_ARRIDX_)+dd(i,j,ci)
         end do
         cc(_ARRIDX_)=(cc(_ARRIDX_)+0.5*dt*ppsum(_ARRIDX_))/(1.+0.5*dt*ddsum(_ARRIDX_)/cc1(_ARRIDX_))
      end do
   end do

   return
   end subroutine patankar_runge_kutta_2
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Fourth-order Patankar-Runge-Kutta scheme
!
! !INTERFACE:
   subroutine patankar_runge_kutta_4(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! This subroutine should become the fourth-order Patankar Runge-Kutta
! scheme, but it does not yet work.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: ppsum(_ARRSZ_),ddsum(_ARRSZ_)
  REALTYPE :: ppsum1(_ARRSZ_),ddsum1(_ARRSZ_)
  REALTYPE :: ppsum2(_ARRSZ_),ddsum2(_ARRSZ_)
  REALTYPE :: ppsum3(_ARRSZ_),ddsum3(_ARRSZ_)
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_
   dd=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do ci=1,nlev
      do i=1,numc
         ppsum(_ARRIDX_)=_ZERO_
         ddsum(_ARRIDX_)=_ZERO_
         do j=1,numc
            ppsum(_ARRIDX_)=ppsum(_ARRIDX_)+pp(i,j,ci)
            ddsum(_ARRIDX_)=ddsum(_ARRIDX_)+dd(i,j,ci)
         end do
         cc1(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum(_ARRIDX_))/(1.+dt*ddsum(_ARRIDX_)/cc(_ARRIDX_))
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,pp,dd)

   do ci=1,nlev
      do i=1,numc
         ppsum1(_ARRIDX_)=_ZERO_
         ddsum1(_ARRIDX_)=_ZERO_
         do j=1,numc
            ppsum1(_ARRIDX_)=ppsum1(_ARRIDX_)+pp(i,j,ci)
            ddsum1(_ARRIDX_)=ddsum1(_ARRIDX_)+dd(i,j,ci)
         end do
         cc1(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum1(_ARRIDX_))/(1.+dt*ddsum1(_ARRIDX_)/cc1(_ARRIDX_))
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,pp,dd)

   do ci=1,nlev
      do i=1,numc
         ppsum2(_ARRIDX_)=_ZERO_
         ddsum2(_ARRIDX_)=_ZERO_
         do j=1,numc
            ppsum2(_ARRIDX_)=ppsum2(_ARRIDX_)+pp(i,j,ci)
            ddsum2(_ARRIDX_)=ddsum2(_ARRIDX_)+dd(i,j,ci)
         end do
         cc1(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum2(_ARRIDX_))/(1.+dt*ddsum2(_ARRIDX_)/cc1(_ARRIDX_))
      end do
   end do

   call right_hand_side(first,numc,nlev,cc1,pp,dd)

   do ci=1,nlev
      do i=1,numc
         ppsum3(_ARRIDX_)=_ZERO_
         ddsum3(_ARRIDX_)=_ZERO_
         do j=1,numc
            ppsum3(_ARRIDX_)=ppsum3(_ARRIDX_)+pp(i,j,ci)
            ddsum3(_ARRIDX_)=ddsum3(_ARRIDX_)+dd(i,j,ci)
         end do
         ppsum(_ARRIDX_)=1./3.*(0.5*ppsum(_ARRIDX_)+ppsum1(_ARRIDX_)+ppsum2(_ARRIDX_)+0.5*ppsum3(_ARRIDX_))
         ddsum(_ARRIDX_)=1./3.*(0.5*ddsum(_ARRIDX_)+ddsum1(_ARRIDX_)+ddsum2(_ARRIDX_)+0.5*ddsum3(_ARRIDX_))
         cc(_ARRIDX_)=(cc(_ARRIDX_)+dt*ppsum(_ARRIDX_))/(1.+dt*ddsum(_ARRIDX_)/cc1(_ARRIDX_))
      end do
   end do

   return
   end subroutine patankar_runge_kutta_4
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: First-order Modified Patankar scheme
!
! !INTERFACE:
   subroutine modified_patankar(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the first-order Modified Patankar-Euler scheme (MPE1) scheme is coded,
! with one evaluation of the right hand side per time step:
! \begin{equation}\label{eq:am:MP}
! \begin{array}{rcl}
!   c_i^{n+1} &=&
! \displaystyle
!  c_i^n
!               + \Delta t \left\{ \sum\limits_{\stackrel{j=1}{j \not= i}}^I
!                p_{i,j}\left(\underline{c}^n\right) \dfrac{c_j^{n+1}}{c_j^n}
!                                         + p_{i,i}\left(\underline{c}^n\right)
!                           - \sum_{j=1}^I d_{i,j}\left(\underline{c}^n\right)
!                                         \dfrac{c_i^{n+1}}{c_i^n} \right\}.
! \end{array}
! \end{equation}
!
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  REALTYPE :: a(1:numc,1:numc),r(1:numc)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_
   dd=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp(i,j,ci)/cc(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc(:,ci))
   end do

   return
   end subroutine modified_patankar
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Second-order Modified Patankar-Runge-Kutta scheme
!
! !INTERFACE:
   subroutine modified_patankar_2(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! Here, the second-order Modified Patankar-Runge-Kutta (MPRK2) scheme is coded,
! with two evaluations of the right hand sides per time step:
!
! \begin{equation}\label{eq:am:MPRK}
!   \left. \begin{array}{rcl}
!     c_i^{(1)} &=&
! \displaystyle
! c_i^n  + \Delta t
! \left\{
! \sum\limits_{\stackrel{j=1}{j \not= i}}^I p_{i,j}\left(\underline{c}^n\right)
! \dfrac{c_j^{(1)}}{c_j^n}
! + p_{i,i}\left(\underline{c}^n\right)
! - \sum_{j=1}^I d_{i,j}\left(\underline{c}^n\right)
! \dfrac{c_i^{(1)}}{c_i^n}
! \right\},
! \\ \\
! c_i^{n+1} &=&
! \displaystyle
! c_i^n  + \dfrac{\Delta t}{2}
!                   \left\{
!                     \sum\limits_{\stackrel{j=1}{j \not= i}}^I
!                       \left(p_{i,j}\left(\underline{c}^n\right)
!                           + p_{i,j}\left(\underline{c}^{(1)}\right)
!                       \right) \dfrac{c_j^{n+1}}{c_j^{(1)}}
!                       + p_{i,i}\left(\underline{c}^n\right)
!                       + p_{i,i}\left(\underline{c}^{(1)}\right)
! \right.\\ \\
!               & &
! \displaystyle
! \left.\phantom{c_i^n  + \dfrac{\Delta t}{2} }
!                   - \sum_{j=1}^I
!                       \left(d_{i,j}\left(\underline{c}^n\right)
!                           + d_{i,j}\left(\underline{c}^{(1)}\right)
!                       \right) \dfrac{c_i^{n+1}}{c_i^{(1)}}
!                   \right\}.
!   \end{array}
!   \right\}
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  REALTYPE :: pp1(1:numc,_ARRSZ_),dd1(1:numc,_ARRSZ_)
  REALTYPE :: a(1:numc,1:numc),r(1:numc)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_ ; pp1=_ZERO_
   dd=_ZERO_ ; dd1=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp(i,j,ci)/cc(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc1(:,ci))
   end do

   call right_hand_side(first,numc,nlev,cc1,pp1,dd1)

   pp=0.5*(pp+pp1)
   dd=0.5*(dd+dd1)

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp(i,j,ci)/cc1(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc1(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc(:,ci))
   end do

   return
   end subroutine modified_patankar_2
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Fourth-order Modified Patankar-Runge-Kutta scheme
!
! !INTERFACE:
   subroutine modified_patankar_4(dt,numc,nlev,cc,right_hand_side)
!
! !DESCRIPTION:
! This subroutine should become the fourth-order Modified Patankar Runge-Kutta
! scheme, but it does not yet work.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,pp,dd)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: pp(:,:,:)
         REALTYPE, intent(out)                :: dd(:,:,:)
      end subroutine right_hand_side
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: pp(1:numc,_ARRSZ_),dd(1:numc,_ARRSZ_)
  REALTYPE :: pp1(1:numc,_ARRSZ_),dd1(1:numc,_ARRSZ_)
  REALTYPE :: pp2(1:numc,_ARRSZ_),dd2(1:numc,_ARRSZ_)
  REALTYPE :: pp3(1:numc,_ARRSZ_),dd3(1:numc,_ARRSZ_)
  REALTYPE :: a(1:numc,1:numc),r(1:numc)
  REALTYPE :: cc1(_ARRSZ_)
  INTEGER  :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!  absolutely essential since not all elements are calculated
   pp=_ZERO_ ; pp1=_ZERO_ ; pp2=_ZERO_ ; pp3=_ZERO_
   dd=_ZERO_ ; dd1=_ZERO_ ; dd2=_ZERO_ ; dd3=_ZERO_

   first=.true.
   call right_hand_side(first,numc,nlev,cc,pp,dd)
   first=.false.

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp(i,j,ci)/cc(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc1(:,ci))
   end do

   call right_hand_side(first,numc,nlev,cc1,pp1,dd1)

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd1(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp1(i,j,ci)/cc1(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc1(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp1(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc1(:,ci))
   end do

   call right_hand_side(first,numc,nlev,cc1,pp2,dd2)

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd2(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp2(i,j,ci)/cc1(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc1(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp2(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc1(:,ci))
   end do

   call right_hand_side(first,numc,nlev,cc1,pp3,dd3)

   pp=1./3.*(0.5*pp+pp1+pp2+0.5*pp3)
   dd=1./3.*(0.5*dd+dd1+dd2+0.5*dd3)

   do ci=1,nlev
      do i=1,numc
         a(i,i)=_ZERO_
         do j=1,numc
            a(i,i)=a(i,i)+dd(i,j,ci)
            if (i.ne.j) a(i,j)=-dt*pp(i,j,ci)/cc1(j,ci)
         end do
         a(i,i)=dt*a(i,i)/cc1(_ARRIDX_)
         a(i,i)=1.+a(i,i)
         r(i)=cc(_ARRIDX_)+dt*pp(i,_ARRIDX_)
      end do
      call matrix(numc,a,r,cc(:,ci))
   end do

   return
   end subroutine modified_patankar_4
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: First-order Extended Modified Patankar scheme
!
! !INTERFACE:
   subroutine emp_1(dt,numc,nlev,cc,right_hand_side)
!
#if 0
! !DESCRIPTION:
! Here, the first-order Extended Modified Patankar scheme for
! biogeochemical models is coded, with one evaluation of the right-hand
! side per time step:
!
! \begin{eqnarray}
!    \lefteqn{\vec{c}^{n + 1} = \vec{c}^n + \Delta t \: \vec{f}(t^n,\vec{c}^n)\prod\limits_{j \in J^n} {\frac{c_j^{n + 1} }{c_j^n}}} \nonumber \\
!    & & \mbox{with } J^n = \left\{ {i:f_i (t^n,\vec{c}^n) < 0,i \in \{1,...,I\}} \right\}
! \end{eqnarray}
!
! This system of non-linear implicit equations is solved in auxiliary subroutine
! findp\_bisection, using the fact this system can be reduced to a polynomial in one
! unknown, and additionally using the restrictions imposed by the requirement of positivity.
! For more details, see \cite{Bruggemanetal2005}.
!
#endif
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)             :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  REALTYPE :: derivative(_ARRSZ_)
  INTEGER  :: ci
  REALTYPE :: pi
!EOP
!-----------------------------------------------------------------------
!BOC
   first=.true.
   call right_hand_side(first,numc,nlev,cc,derivative)
   first=.false.

   do ci=1,nlev
      call findp_bisection(numc, cc(:,ci), derivative(:,ci), dt, _POINT_NEG_NINE_, pi)
      cc(:,ci) = cc(:,ci) + dt*derivative(:,ci)*pi
   end do

   return
   end subroutine emp_1
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Second-order Extended Modified Patankar scheme
!
! !INTERFACE:
   subroutine emp_2(dt,numc,nlev,cc,right_hand_side)
!
#if 0
! !DESCRIPTION:
! Here, the second-order Extended Modified Patankar scheme for
! biogeochemical models is coded, with two evaluations of the right-hand
! side per time step:
!
! \begin{eqnarray}
!   \vec{c}^{(1)} & = & \vec{c}^n  + \Delta t \: \vec{f}(t^n ,\vec{c}^n )\prod\limits_{j \in J^n } {\frac{{c_j^{(1)} }}{{c_j^n }}} \nonumber \\
!   \vec{c}^{n + 1} & = & \vec{c}^n  + \frac{{\Delta t}}{2}\left( {\vec{f}(t^n ,\vec{c}^n ) + \vec{f}(t^{n + 1} ,\vec{c}^{(1)} )} \right)\prod\limits_{k \in K^n } {\frac{{c_k^{n + 1} }}{{c_k^{(1)} }}}
! \end{eqnarray}
!
! where
!
! \begin{eqnarray}
!   J^n & = & \left\{ {i:f_i (t^n ,\vec{c}^n ) < 0, i \in \{ 1,...,I\} } \right\} \nonumber \\
!   K^n & = & \left\{ {i:f_i (t^n ,\vec{c}^n ) + f_i (t^{n+1} ,\vec{c}^{(1)} ) < 0, i \in \{ 1,...,I\} } \right\}.
! \end{eqnarray}
!
! The first step is identical to a step with the first-order EMP scheme. The second step mathmatically identical to
! a step with the first-order scheme if we rewrite it as
!
! \begin{eqnarray}
!   \lefteqn{\vec{c}^{n + 1} = \vec{c}^n + \Delta t \: \vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)})\prod\limits_{k \in
!     K^n} {\frac{c_k^{n + 1} }{c_k^n }}} \nonumber \\
!   & & \mbox{with }\vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)}) = \frac{1}{2}\left( {\vec{f}(t^n,\vec{c}^n) + \vec{f}(t^{n + 1},\vec{c}^{(1)})} \right)\prod\limits_{k \in K^n} {\frac{c_k^n }{c_k^{(1)} }}.
! \end{eqnarray}
!
! Therefore, this scheme can be implemented as two consecutive steps with the first-order scheme, the second using
! $\vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)})$. The non-linear problem of each consecutive step is solved
! in auxiliary subroutine findp\_bisection.
!
! For more details, see \cite{Bruggemanetal2005}.
!
#endif
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: dt
   INTEGER, intent(in)                 :: numc,nlev
!
! !INPUT/OUTPUT PARAMETER:
   REALTYPE, intent(inout)              :: cc(:,:)

   interface
      subroutine right_hand_side(first,numc,nlev,cc,rhs)
         USE ISO_C_BINDING
         LOGICAL, intent(in)                  :: first
         INTEGER, intent(in)                  :: numc,nlev
         REALTYPE, intent(in)                 :: cc(:,:)
         REALTYPE, intent(out)                :: rhs(:,:)
      end
   end interface
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
! !LOCAL VARIABLES:
  LOGICAL  :: first
  INTEGER  :: _ARRIDX_
  REALTYPE :: pi, rhs(_ARRSZ_), cc_med(_ARRSZ_), rhs_med(_ARRSZ_)
!EOP
!-----------------------------------------------------------------------
!BOC
   first=.true.
   call right_hand_side(first,numc,nlev,cc,rhs)
   first=.false.

   do ci=1,nlev
      call findp_bisection(numc, cc(:,ci), rhs(:,ci), dt, _POINT_NEG_NINE_, pi)
      cc_med(:,ci) = cc(:,ci) + dt*rhs(:,ci)*pi
   end do

   call right_hand_side(first,numc,nlev,cc_med,rhs_med)

   do ci=1,nlev
      rhs(:,ci) = 0.5 * (rhs(:,ci) + rhs_med(:,ci))

      ! Correct for the state variables that will be included in 'p'.
      do i=1,numc
         if (rhs(_ARRIDX_) .lt. 0.) rhs(:,ci) = rhs(:,ci) * cc(_ARRIDX_)/cc_med(_ARRIDX_)
      end do

      call findp_bisection(numc, cc(:,ci), rhs(:,ci), dt, _POINT_NEG_NINE_, pi)

      cc(:,ci) = cc(:,ci) + dt*rhs(:,ci)*pi
   end do ! ci (z-levels)

   return
   end subroutine emp_2
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Calculation of the EMP product term 'p'
!
! !INTERFACE:
   subroutine findp_bisection(numc, cc, derivative, dt, accuracy, pi)
!
#if 0
! !DESCRIPTION:
! Auxiliary subroutine for finding the Extended Modified Patankar
! product term $p$ with the bisection technique.
!
! This subroutine solves the non-linear problem
!
! \begin{eqnarray}
!    \lefteqn{\vec{c}^{n + 1} = \vec{c}^n + \Delta t \: \vec{f}(t^n,\vec{c}^n)\prod\limits_{j \in J^n} {\frac{c_j^{n + 1} }{c_j^n}}} \nonumber \\
!    & & \mbox{with } J^n = \left\{ {i:f_i (t^n,\vec{c}^n) < 0,i \in \{1,...,I\}} \right\}
! \end{eqnarray}
!
! using the fact that it can be reduced to the problem of finding the root of a polynomial
! in one unknown $p: = \prod\limits_{j \in J^n}{{c_j^{n + 1}}/{c_j^n}}$:
!
! \begin{equation}
!   g(p) = \prod\limits_{j \in J^n} {\left( {1 + \frac{\Delta t \: f_j(t^n,\vec{c}^n)}{c_j^n }p} \right)} - p = 0,
! \end{equation}
!
! with
!
! \begin{equation}
!   J^n = \left\{ {i:f_i (t^n,\vec{c}^n) < 0,i \in \{1,...,I\}} \right\},
! \end{equation}
!
! Additionally, it makes use of the the positivity requirement $\vec{c}^{n+1}_i>0\ \forall\ i$, which
! imposes restriction
!
! \begin{equation}
!   p \in \left( 0, \min \left( {1,\mathop {\min }\limits_{j \in J^n} \left( { -
!                   \frac{c_j^n }{\Delta t \: f_j (t^n,\vec{c}^n)}} \right)} \right) \right).
! \end{equation}
!
! It has been proved that there exists exactly one $p$ for which the above is
! true, see \cite{Bruggemanetal2005}.
! The resulting problem is solved using the bisection scheme, which is guaranteed to converge.
!
#endif
!
! !USES:
   implicit none
!
! !INPUT PARAMETERS:
   INTEGER, intent(in)   :: numc
   REALTYPE, intent(in)  :: cc(:), derivative(:)
   REALTYPE, intent(in)  :: dt, accuracy
!
! !OUTPUT PARAMETER:
   REALTYPE, intent(out) :: pi
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
! !LOCAL VARIABLES:
   REALTYPE :: pileft, piright, fnow
   REALTYPE :: relderivative(1:numc)
   INTEGER  :: iter, i, potnegcount
!EOP
!-----------------------------------------------------------------------
!BOC
! Sort the supplied derivatives (find out which are negative).
   potnegcount = 0
   piright = 1.
   do i=1,numc

      if (derivative(i).lt.0.) then
!        State variable could become zero or less; include it in the
!        J set of the EMP scheme.
         if (cc(i).eq.0.) write (*,*) "Error: state variable ",i," is zero and has negative derivative!"
         potnegcount = potnegcount+1
         relderivative(potnegcount) = dt*derivative(i)/cc(i)

!        Derivative is negative, and therefore places an upper bound on pi.
         if (-1./relderivative(potnegcount).lt.piright) piright = -1./relderivative(potnegcount)
     end if

   end do

   if (potnegcount.eq.0) then
!     All derivatives are positive, just do Euler.
      pi = 1.0
      return
   end if

   pileft = 0.      ! polynomial(0) = 1

!  Determine maximum number of bisection iterations from
!  requested accuracy.
!  maxiter = -int(ceiling(dlog10(accuracy)/dlog10(2.D0)))

   do iter=1,20
!     New pi to test is middle of current pi-domain.
      pi = 0.5*(piright+pileft)

!     Calculate polynomial value.
      fnow = 1.
      do i=1,potnegcount
         fnow = fnow*(1.+relderivative(i)*pi)
      end do

      if (fnow>pi) then
!        Polynomial(pi)>0; we have a new left bound for pi.
         pileft = pi
      elseif (fnow<pi) then
!       Polynomial(pi)<0; we have a new right bound for pi.
        piright = pi
      else
!       Freak occurrence: polynomial(pi)=0, we happened to pinpoint
!       the exact pi.
        exit
      end if
!     Check if we now pi accurately enough (accuracy refers to the
!     number of decimals we know).
      if ((piright-pileft)/pi<accuracy) exit
   end do

!  Low pi values imply very large negative relative derivative. This happens
!  for stiff systems (or very high delta_t), and for non-positive systems.
!  Then EMP is not suitable (it will stall state variable values), so warn user.
   if (pi.lt.1.d-4) then
     write (*,*) "Warning: small pi=",pi," in Extended Modified Patankar slows down system!"
!    write (*,*) "relative derivatives: ",derivative(:)*dt/cc(:)
!    write (*,*) "You system may be stiff or non-positive, or you time step is too large."
!    stop "ode_solvers::findp_bisection"
   end if

   return

   end subroutine findp_bisection
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Matrix solver
!
! !INTERFACE:
   subroutine matrix(n,a,r,c)
!
! !DESCRIPTION:
! This is a Gaussian solver for multi-dimensional linear equations.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   INTEGER, intent(in)                 :: n
!
! INPUT/OUTPUT PARAMETERS:
  REALTYPE                             :: a(:,:),r(:)
!
! OUTPUT PARAMETERS:
  REALTYPE, intent(out)                :: c(:)
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  INTEGER  :: i,j,k
!EOP
!-----------------------------------------------------------------------
!BOC
   do i=1,n
      r(i)=r(i)/a(i,i)
      do j=n,i,-1
         a(i,j)=a(i,j)/a(i,i)
      end do
      do k=i+1,n
         r(k)=r(k)-a(k,i)*r(i)
         do j=i+1,n
            a(k,j)=a(k,j)-a(k,i)*a(i,j)
         end do
      end do
   end do

   do i=n,1,-1
      c(i)=r(i)
      do j=i+1,n
         c(i)=c(i)-a(i,j)*c(j)
      end do
   end do

   return
   end subroutine matrix

END MODULE
!EOC

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------
